Slightly Advanced Introduction to Mathematica.

Below is a slightly advanced introduction to some concepts in Mathematica.  You may be familiar with most of them.

Simple Things:

Note that there is a difference between the two types of answers below.  Mathematica tries to handle things symbolically if it can.  Often, we don't want that.

Sin[π/3]
Sin[π/3.]

Integrals:

Mathematica can do integrals symbolically.   Here's an integral you probably don't want to do by hand:

Integrate[Sin[x]/x,{x,0,Infinity}]

Mathematica can do this integral numerically.  But note that  Mathematica fails miserably, unless we give it some help.

NIntegrate[Sin[x]/x,{x,0,Infinity}]

NIntegrate[Sin[x]/x,{x,0,Infinity},Method-> Oscillatory]

The point is that Mathematica can also do integrals that have no hope of being done analytically:

NIntegrate[((BesselJ[0,x])^2) Sqrt[(1/(1+(x^2)))] ,{x,0,1}]

Solving Systems of Equations

Sometimes you might like to solve a system of equations analytically.   Consider the standard problem of scattering an electron of energy ϵ  from a square barrier of height v.   You will recall that we have wavefunctions in each of the regions:

ψ1 = Exp[I k x ] + a Exp[-I k x ]
ψ2 =b Exp[I κ x]   + c  Exp[- I κ x]
ψ3 = d Exp[I k x]

We need their derivatives :

dψ1 = D[ψ1,x];
dψ2 = D[ψ2,x];
dψ3 = D[ψ3,x];

Now we need to evaluate them at x=0 and x=L.

ψ1at0=ψ1/. x-> 0;
ψ2at0=ψ2/. x-> 0;

ψ2atL=ψ2/. x-> L;
ψ3atL=ψ3/.x-> L;

dψ1at0=D[ψ1,x]/. x-> 0;
dψ2at0=D[ψ2,x]/. x-> 0;

dψ2atL=D[ψ2,x]/. x-> L;
dψ3atL=D[ψ3,x]/.x-> L;

We have all of the pieces that we need.   Now we can put them into a system of equations.  Note the use of the "==" in systems of equations.

eqns={ψ1at0==ψ2at0, dψ1at0==dψ2at0, ψ2atL==ψ3atL, dψ2atL==dψ3atL}

Solve the system:

solution=Solve[eqns,{a,b,c,d}]

t=(d/.solution)

t=(d/.solution)[[1]]

v0= 32.;
k0=(Sqrt[ϵ]) ;
κ0= (Sqrt[ϵ - v0]);
tplot=t/. {k-> k0, κ-> κ0, L-> 1.}

Plot[((Abs[tplot])^2),{ϵ, .1,100}]

ψ1sol=(ψ1/. solution)[[1]] ;
ψ2sol=(ψ2/. solution)[[1]];
ψ3sol=(ψ3/. solution)[[1]];
ψsol=Piecewise[{{ ((Abs[ψ1sol])^2),x< 0},{ ((Abs[ψ2sol])^2),x>0 && x<L},{ ((Abs[ψ3sol])^2),x>L}}];

v0= 32.;
ϵ0=55.;
k0=(Sqrt[ϵ0]) ;
κ0= (Sqrt[ϵ0 - v0]);
ψplot=ψsol/. {k-> k0, κ-> κ0, L-> 1.}

Plot[ψplot,{x,-4,4}]

Vectors and Matrices:

The fundamental item in Mathematica is the "list."   This is much like a vector. It consists of a set of elements  enclosed in curly brackets:  "{a1, a2, a3...}".    We can construct a list either by hand, or by having Mathematica perform some action repeatedly using the Table command.  The latter is much like a "do-loop" in Fortran or a "for-next" loop in Basic.

vector1={1,2,3,2}
vector2=Table[n!,{n,1,20}]

When working with vectors and matrices, it is important to know when to tell Mathematica to shut up.  If you put a semicolon on the end of a line, then the calculation is performed, but the output is suppressed.

vector1={1,2,3,4};
vector2=Table[n!,{n,1,20}];

If you want a pretty version of the vector, then use MatrixForm.

MatrixForm[vector2]

To refer to an element in the list, we just double square brackets.  For example, the 10th element in vector2 is vector2[[10]]:

vector2[[10]]

A matrix is just a list of lists.  Again, we can write it out by hand or have Mathematica create it using a Table command.

matrix1={{1,2,3,4},{5,0,7,8},{5,4,0,2},{7,6,5,4}};
matrix2= Table[  Sin[i*j],{i,1,4},{j,1,4}];
MatrixForm[matrix1]
MatrixForm[matrix2]

Again, we can address an individual element using the double square brackets, but now we have to have two indices:

matrix1[[2,3]]

Inverse[matrix1]

Mathematica will find exact solutions, often when we don't really want them.  Please note the difference between the two commands below and in their output.  Students often write slow Mathematica code because they forget to tell Mathematica to work numerically and not symbolically.

Inverse[matrix2]
Inverse[N[matrix2]]

Mathematica can solve linear sets of equations, and quickly find eigenvalues and eigenvectors.

LinearSolve[matrix1,vector1]

matrix3=N[matrix1];
Eigenvalues[matrix3]
Eigenvectors[matrix3]

Are the eigenvectors normalized?  Let's check:

eigvec=Eigenvectors[matrix3];
Table[eigvec[[i]].Conjugate[eigvec[[i]] ],{i,1,4}]

Pretty darn close!  Let's chop off anything less that (10^(-10)).

Table[Chop[eigvec[[i]].Conjugate[eigvec[[i]] ]],{i,1,4}]

Combining Plots

The simplest way to combine 2D plots is to give them as a list to the plot command.

Plot[{x,(x^2),(x^3),(x^4)},{x,0,1}]

If you name individual plots, you can also combine them using Show.

plt1=Plot[x,{x,0,1}]
plt2=Plot[(x^2),{x,0,1}]
Show[plt1,plt2]

plt3=Plot3D[Sin[x y],{x,-π,π},{y,-π,π},PlotPoints->40]
plt4=Plot3D[Cos[x y],{x,-π,π},{y,-π,π},PlotPoints->40]
Show[plt3,plt4]